/***************************************************************************
 *   Copyright (C) 2007 by Alexandre Masselot,,,   *
 *   alex@alex-laptop   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU Library General Public License as       *
 *   published by the Free Software Foundation; either version 2 of the    *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this program; if not, write to the                 *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#ifndef BOLMASSINDEX_H
#define BOLMASSINDEX_H

#include "deltamass.h"
#include "massindexanchor.h"
#include <vector>
namespace bol {

/**
  An index to anchor register mass + fasta acess to data within a tolerance

  This template class is instantiated for example like MassIndex<MassIndexAnchorDummy, vector<const MassIndexAnchorDummy*>::const_iterator> where MassIndexAnchorDummy has a mass() method

  THe idea is to have a sorted vector, a list of *MassIndexAnchor (typically a mass + other info). For a faster access, we can also set a vector of iterators to position into this list  with a predefined bin

  Warning: after the list is filled, you must call build_index to sort the vector and build the index


  @author Alexandre Masselot,,, <alex@alex-laptop>
 */
  template <class TAnchor, class TAnchorIterator> class MassIndex{
    private:
      vector<const TAnchor*> m_list;
      vector<TAnchorIterator> m_index;
      double m_min, m_max, m_step;
      int m_nbbins;
      const DeltaMass *m_tolerance;

      ///@todo I don't like all the rounding error that can occur here...
      int m2i(const double m)const{
        if(m<=m_min) return 0;
        if(m>=m_max) return nbbins();
        return (int)((m-m_min)/m_step);
      }
      double i2m(const int i)const{
        if(i<=0) return m_min;
        if(i>nbbins()) return m_max;
        return m_min + i*m_step;
      }
      static bool mianchorLT(const TAnchor* m1, const TAnchor* m2){
        return m1->mass()<m2->mass();
      }

    public:

    /**
     * build a index with a given number of bins
      
      This template class is instantiated for example like MassIndex<MassIndexAnchorDummy, vector<const MassIndexAnchorDummy*>::const_iterator> where MassIndexAnchorDummy has a mass() method

     * @param n the number of bins
     */
      MassIndex<TAnchor, TAnchorIterator>(const unsigned int n){
        m_nbbins=n;
        m_min=m_max=0;
      }

      ~MassIndex(){}

      /**
       * add an anchor. class TAnchor must provide a mass() parameter;
       * @param anchor
       */
      void addAnchor(const TAnchor* anchor){
        m_list.push_back(anchor);
      }

      /**
       * sort the registrated anchors and build the index
       */
      void buildIndex(){
        sort(m_list.begin(), m_list.end(), mianchorLT);
        m_index.assign(nbbins()+1, m_list.begin());
        m_index[nbbins()]=m_list.end();
        if(m_list.size()==0){
          for (int i=0; i<nbbins(); i++){
            m_index[i]=m_list.begin();
          }
          return;
        }
        m_min=m_list[0]->mass()-.1;
        m_max=m_list[m_list.size()-1]->mass()+.1;
        m_step=(m_max-m_min)/m_nbbins;
        TAnchorIterator it=m_list.begin();
        double itmass=m_min;
        int iti=0;
        for(iti=0; iti<m2i((*m_list.begin())->mass()) && iti<=nbbins(); iti++){
          m_index[iti]=m_list.begin();
        }
        for(; it<m_list.end() && iti<m_nbbins; it++){
          if((*it)->mass()>=itmass){
            m_index[iti]=it;
            iti++;
            itmass=i2m(iti);
            if(it>m_list.begin())it--;
          }
        }
        for(;iti<m_nbbins;iti++){
          m_index[iti]=m_list.end();
        }
      }
      /**
       * browse the index and set iterators  with the elements close enough to m (with tolerance())
       * @param mass the original mass
       * @param first it greater or equal to first
       * @param last  and strictly lower than last fullfil the tolerance
       * @return the number of matches (vanchors.size())
       */
      int getAnchors(const double mass,
                     TAnchorIterator & first,
                     TAnchorIterator & last
                    )const
      {
        double inf, sup;
        tolerance()->range(mass, inf, sup);

        if(inf<m_min){
          //it out of the index, must scan from begin
          first=m_list.begin();
        }else{
          int ifirst=m2i(inf);
          first=m_index[ifirst];
        }
        if(sup>m_max){
          last=m_list.end();
        }else{
          int ilast=m2i(sup)+1;
          if(ilast>=m_index.size()){
            last=m_list.end();
          }else{
            last=m_index[ilast];
          }
        }
        for(; first<m_list.end() && (*first)->mass()<inf; ++first);
        for(; last>m_list.begin() && (*(last-1))->mass()>sup; --last);
        return last-first;
      }

      /**
       * iterators to the whole mass list
       */
      TAnchorIterator begin() {
        return m_list.begin();
      }
      /**
       * iterators to the whole mass list
       */
      TAnchorIterator end() {
        return m_list.end();
      }

      /**
       * set iterators to begin and end of the whole mass list
       * @param begin is the first mass
       * @param end  is the end of the mass list
       * @return the total number of masses
       * @note this total may be more than the number spectra if some spectra have several possible charge state
       */
      int getAnchors(TAnchorIterator & begin,
                     TAnchorIterator & end
                    )const
      {
	begin = m_list.begin();
	end = m_list.end();
        return end-begin;
      }

      /**
       *
       * @param theValue  the tolerance (used to getAnchors)
       */
      void setTolerance(const DeltaMass* theValue) {
        m_tolerance = theValue;
      }
  

      /**
       *
       * @return the tolerance (used to getAnchors)
       */
      const DeltaMass* tolerance() const {
        return m_tolerance;
      }
  
      /**
       * output the list of filled index + anchjor object
       */
      template <class U, class V> friend ostream& operator<<(ostream& out, const MassIndex<U, V>& idx){
        out<<"list"<<endl;
        for(V it=idx.m_list.begin(); it<idx.m_list.end(); it++){
          cout<< **it<<endl;
        }

        out<<"index ["<<idx.m_min<<":"<<idx.m_max<<"]/"<<idx.m_step<<endl;
    //vector<const MassIndexAnchor*>::const_iterator last=*(idx.m_index.begin()-1);
        for (int i=0; i<idx.m_index.size(); i++){
      //if(idx.m_index[i]!=last){
          out<<i<<"("<<idx.i2m(i)<<")"<<"\t";

          if(idx.m_index[i]==idx.m_list.end()){
            out<<"END";
          }else{
            out<< **idx.m_index[i];
          }
          out<<endl;
/*      }
          last=idx.m_index[i];*/
        }
        return out;
      }

      int nbbins() const {
        return m_nbbins;
      }
  
  };
}

#endif
